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CN 1 Introduction 

r I This file is the Sweave documentation for the examples provided in |Flegal and| 
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1.1 Batch Means 



The following function is required to calculate the variance estimators using 
batch means (BM) and overlapping batch means (OLBM). We are thankful to 
Murali Haran whom wrote the original function to implement BM which we 
have expanded here. 

> id <- function (x) return (x) 



> mcse <- functionCvals , bs = "sqroot" , g = id, meth = "BM" , warn = FALSE) { 



+ N <- length(vals) 



+ if (N < 1000) { 
+ if (warn) 



+ cat("WARNING: too few samples (less than 1000)\n") 

+ if (N < 10) 

+ return(NA) 

^ + if (bs == "sqroot") { 

+ b <- floor (sqrt(N)) 

+ a <- floor (N/b) 

+ } 

+ else if (bs == "cuberoot") { 

+ b <- floor (N" (1/3)) 

+ a <- floor (N/b) 

+ } 

+ else { 

+ stopifnot(is.numeric(bs)) 
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+ 




bmse <- sqrt (sigmsLhatsq/N) 


+ 




return (bmse) 




} 






if 


(meth == "OBM") { 


+ 




a <- N - b + 1 


+ 




Ys <- sapply(l:a, function(k) return (mean (g (vals [k: (k + 


+ 




b - 1)])))) 


+ 




mubat <- meaii(Ys) 


+ 




sigmahatsq <- N * b * sum((Ys - muhat) ~2) / (a - l)/a 


+ 




bmse <- sqrt (sigmahatsq/N) 


+ 




return(bmse) 


+ 


} 






else { 






stopC'method specified invalid (meth=", meth, ")") 




} 




+ } 







2 Normal AR(1) Markov Chains 

Consider the normal AR(1) time series defined by 

Xn+l = pXn + en (1) 

where the e„ are i.i.d. N(0,1) and p < 1. This Markov chain has invariant 
distribution N (0, 1/(1 - p^)). 

In our example, we perfomred calulcations for p e {0.5,0.95}. 

2.1 Markov Chain Sampler 

The following chunk of code gives general functions needed to sample from 
0. Here we have a function that provides an observation from the chain, an 
observation q steps ahead with a defualt of one step ahead, and p observations 
from the chain. 

> arl <- function(m, rho , tau) { 
+ rho * m + rnormd, 0, tau) 
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+ } 

> arl.q <- function ("in, rho, tau, q = 1) { 
+ for a in l:q) { 

+ m <- rho * m + rnormCl, 0, tau) 

+ } 
+ m 

+ } 

> arl.gen <- function(mc , p, rho, tau, q = 1) { 
+ loc <- leiigth(mc) 

+ junk <- double (p) 

+ mc <- append (mc, junk) 

+ for a in l:p) { 

+ j <- i + loc - 1 

+ mc[(j + 1)2 <- arl(mc[j] , rho, tau) 

+ } 

+ return (mc) 
+ } 

2.2 Additional Functions 

The following are additional functions necessary for later calculations. The first 
calculates the estimated first and third quartile while the second calculates the 
associated MCSE via subsampling. Notice this function is similar to the function 
necessary for OLBM. 



> quant <- function ("inputs { 

+ quant He (input, prob = c(0.25 , 0.75), type = 1) 
+ } 

> subsampling <- function(vals) { 
+ N <- length(vals) 

+ b <- floor (sqrt(N)) 

+ a <- N - b + 1 

+ Ys <- sapplyd :a, function(k) return(quant(vals [k: (k + b - 
+ 1)]))) 

+ muhat <- apply (Ys, 1, mean) 

+ sigmahatsq <- N * b * apply((Ys - muhat)~2, 1, sum)/(a - 
+ l)/a 

+ bmse <- sqrt(sigmahatsq/N) 

+ return(bmse) 
+ } 



2.3 Simulation Settings and Calculations 

In this next chunk of code, we first give the simulation settings used throughout 

the paper. We then generate the two Markov chains for p € {0.5,0.95} and 
calculate the corresponding estimates and MCSEs. 
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> 22 <- 2000 

> iter <- seq(l, n) 

> crit.bm <- qt (0.9, (sqrt(iter) - 1)) 

> crit.obm <- qt(0.9, (iter - sqrt(iter))) 

> set. seed (1976) 

> rhol <- 0.5 

> rho2 <- 0.95 

> chainl <- arl .gen(l , (n - 1) , rhol, 1) 

> meanl <- cumsum(chainl)/seq(along = chainl) 

> bm.estl <- sapply(l:n, function(k) returnCmcse (chainl [1 :kj , meth = "BM"))) 

> obm.estl <- sapply(l:n, function(k) return(mcse (chainl [1 : k] , 
+ meth = "OBM"))) 

> quartilel <- sapply(l:n, function(k) return(quant (chainl [1 :k] )) ) 

> q.mcsel <- sapply(l:n, function(k) return(subsampling(chainl[l:k]))) 

> chain2 <- arl .gen(l , (n - 1) , rho2, 1) 

> mean2 <- cumsum(chain2)/seq(along = chain2) 

> bm.est2 <- sapply(l:n, function(k) return (mcse (chain2 [1 : k] , meth = "BM"))) 

> obm.est2 <- sapply(l :n, function(k) return(mcse(chain2[l:kj , 

+ meth = "OBM"))) 

> quartile2 <- sapply(l:n, function(k) return(quant(chain2[l:k]))) 

> q.mcse2 <- sapply(l:n, function(k) return(subsampling(chain2[l:k]))) 

2.4 Initial Examination of Output 

The following chunk of code will create plots for the initial examination of 
output. The plots are also repeated in the document. 

> rho = rhol 

> chain <- chainl 

> mean <- meanl 

> par (mf row = c(3, 1) , mar = c(3, 4, 4, 2)) 

> ts. plot (chain, main = "Time-Series vs. Iteration", xlab = "", 

+ ylab = xlim = c(0, n) ) 

> abline(h = 2 * sqrt(l/(l - rho~2))) 

> abline(h = -2 * sqrt(l/(l - rho~2))) 

> acf (chain, main = "Autocorrelation vs. Lag", ylab = "", xlab = "") 

> ts.plot(mean, main = "flu2222i22g' Average vs. Iteration" , xlab = "", 
+ ylab = "", Iwd = 2, xlim = c(0, n)) 

> abline(h = 0) 

> par (mf row = c(l, 1) , mar = c(5, 4, 4, 2)) 
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> rho = rho2 

> chain <- cbain2 

> mean <- mean2 

> par (mf row = c(3, 1) , mar = c(3, 4, 4, 2)) 

> ts .plot (chain, main = "Time-Series vs. Iteration" , xlab = "", 
+ ylah = "", xlim = c(0, n)) 

> ablineCh = 2 * sqrt(l/(l - rho'2))) 

> ablineCh = -2 * sqrt(l/(l - rho'2))) 

> acf (chain, main = "Autocorrelation vs. Lag", ylab = "", xlab = "") 

> ts .plot(mean, main = "Running Average vs. Iteration" , xlab = "", 
+ ylab = "", Iwd = 2, xlim = c(0, n) ) 

> abline(h = 0) 

> par (mf row = c(l, 1), mar = c(5, 4, 4, 2)) 
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Time-Series vs. Iteration 




2.5 Running MCSEs for Expectations 

The following chunk of code creates the plot of the running MCSEs and run- 
ning estimates for the expectations with confidence bounds. Again, the plot is 
contained in the document. 

> rho = rhol 

> chain <- chainl 

> mean <- meanl 

> u.obm <- mean + crit.obm * obm.estl 

> l.obm <- mean - crit.obm * obm.estl 

> u.bm <- mean + crit.bm * bm.estl 

> l.bm <- mean - crit.bm * bm.estl 

> ts .plotCmean, main = "Running Average" , xlab = "Iteration" , ylab = "", 
+ yiiin = c(min(l.obm[10:n]) , max(u.obm[10:n] )) , Iwd = 2, xlim = c(0, 
+ n)) 

> abline(h = 0) 

> points (iter, u.obm, type = "1", Ity = 4, Iwd = 2) 

> pointsdter, l.obm, type = "1", Ity = 4, Iwd = 2) 



6 



Running Average 

o 




Iteration 



> rho = rho2 

> chain <- chain2 

> mean. <- meaji2 

> u.obm <- mean + crit.obm * obm.est2 

> l.obm <- mean - crit.obm * obm.est2 

> u.bm <- mean + crit.bm * bm.est2 

> l.bm <- mean - crit.bm * bm.est2 

> ts . plot (mean , main = "Running Average" , xlab = "Iteration", ylab = "", 
+ ylim = c(min(l.obm[10:n]) , max(u.obm[10:n] )) , Iwd = 2, xlim = c(0, 
+ n)) 

> abline(h = 0) 

> pointsdter, u.obm, type = "1", Ity = 4, Iwd = 2) 

> pointsdter, l.obm, type = "1", Ity = 4, Iwd = 2) 
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Running Average 




2.6 Running QuEirtile Plots 

Here is the code for the running quartile plots with and without confidence 
bounds for p G {0.5,0.95}. There are a total of four chunks of code here, one 
for each plot. 

> rho = rhol 

> chain <- chainl 

> quartiles <- quartilel 

> u.sub <- quartiles + rbind(cr it. obm, crit.obm) * q.mcsel 

> I. sub <- quartiles - rbind(crit .obm, crit.obm) * q.mcsel 

> ts. plot (t (quartiles) , main = "Running Quartiles" , xlab = "Iteration" , 
+ ylab = "", Iwd = 2, ylim = c (max (u . sub [ , 100:n]), min(l.sub[, 

+ lOO-.n])), xlim = c(0, n)) 

> abline(h = qnorm(0.25, 0, sqrt(l/(l - rho'2)))) 

> abline(h = qnorm(0.75, 0, sqrt(l/(l - rho~2)))) 
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> ts. plot (t (quartiles) , main = "Running Quartiles" , xlab = "Iteration" , 



+ 


ylab = "", Iwd = 2, 


ylim = c(max(u. 


sub[. 


100 


■n]). 


min 


+ 


100:n])), xlim 


= c(0, n)) 










> 


points (iter , t(u. sub[l , 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 


> 


points (iter , t(u.sub[2, 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 


> 


points (iter , t(l.sub[l, 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 


> 


points (iter , t(l.sub[2, 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 



> abline(h = qnorm(0.25, 0, sqrt(l/(l - rho'2)))) 

> abline(b = qnorm(0.75, 0, sqrt(l/(l - rbo''2)))) 
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> rho = rho2 

> chain <- chain2 

> quartiles <- quartile2 

> u.sub <- quartiles + rhindCcr it .ohm, crit.obm) * q.mcse2 

> l.sub <- quartiles - rbind (cr it . obm, crit.obm) * q.mcse2 

> ts. plot (t (quartiles) , main = "Running Quartiles" , xlab = "Iteration", 
+ ylab = "", Iwd = 2, ylim = c (max (u . sub [ , 100:n]), min(l.sub[, 

+ 100:n])), xlim = c(0, n)) 

> abline(b = qnorm(0.25, 0, sqrt(l/(l - rbo~2)))) 

> abline(b = qnorm(0.75, 0, sqrt(l/(l - rbo''2)))) 
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> ts. plot (t (quartiles) , main = "Running Quartiles" , xlab = "Iteration" , 



+ 


ylab = "", Iwd = 2, 


ylim = c(max(u. 


sub[. 


100 


■n]). 


min 


+ 


100:n])), xlim 


= c(0, n)) 










> 


points (iter , t(u. sub[l , 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 


> 


points (iter , t(u.sub[2, 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 


> 


points (iter , t(l.sub[l, 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 


> 


points (iter , t(l.sub[2, 


]), type = "1", 


Ity 


= 4, 


Iwd 


= 2) 



> abline(h = qnorm(0.25, 0, sqrt(l/(l - rho'2)))) 

> abline(b = qnorm(0.75, 0, sqrt(l/(l - rbo''2)))) 
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2.7 Calculations for Expectations 

The following are calculations reported in the paper based on n = 2000 iterations 
in the chain. 

> signif (meanl [n] , 6) 
[1] -0.0336436 

> signif (crit.obmEn] * obm.estl[n] , 6) 
[1] 0.0556022 

> signif (meaii2 [n] , 6) 
[1] -0.506755 

> signif (crit.obm[n] * obm.est2[n], 6) 
[1] 0.450564 

> signif (obm . est2 [n] , 6) 
[1] 0.351458 
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Here the implementation of fixed width procedures using p = 0.95. Here we 
have added 198000 iterations to the chain for a total of 200000 to spead up the 
processing time. If one desired a level of precision smaller than 0.1 it would 
probably be necessary to add additional iterations to the chain. 

> chains <- arl.gen(cbaiii2, 198000, rho2, 1) 

> chain <- chain2 

> half <- crit.obm[n] * obm. est2 [n] 

> N <- n 

> while (half + 1/N > 0.1) { 
+ N <- N + 1000 

+ b <- floor (sqrt(N)) 

+ t.OBM <- qt (0.9, (N - b + 1)) 

+ est. OBM <- mcse (chain [1 : N] , metb = "OBM") 

+ half <- t.OBM * est. OBM 

+ } 

> N 

[1] 60000 

> half 

[1] 0.0996075 

> mean ("chain [1 : NJ ) 
[1] -0.04415015 

> signif (mean(chain[l :N] ) , 6) 
[1] -0.0441502 

2.8 Calculations for Quartiles 

Here are the same calculations for the quartiles. First the caculations based on 
n = 2000 iterations for the quartiles. 

> signif (quartilel [, n] , 6) 

257. 757. 
-0.816573 0.777758 

> signif (crit . obm[n] * q.mcsel [, n] , 6) 

257. 757. 
0.0685046 0.0648736 

> signif (quartile2[, n] , 6) 



13 



25% 757. 
-2.73625 1.77970 



> signif(crit.obm[n] * q.mcse2[, n] , 6) 

257. 757. 
0.480563 0.466388 



Then the fixed width calculations for quantilcs with p = 0.95. First for 
individual CIs, the with a Bonferonni correction. These calculations are not 
included in the paper and are commented out since the computational time is 
approximately an hour. The code is included for future reference. 

chain <- chain2 

half <- inax(crit.obm[n]*q.mcse2[,n]) 

N <- n 

while (half + 1 / N > .1){ 
M <- N + 2000 

b <- floor (sqrt(M)) # batch size 
t.sub <- qt(.9, (N - b + 1)) 
est. sub <- subsainpling(chain[l :N] ) 
half <- max(t.sub*est.sub) 

> 
N 

half 

quant (chain [1 : M] ) 
t . sub*est . sub 



##### 

# Then for Bonferonni Correction 
##### 



t.sub <- t.sub <- qt(.9875, (N - b + 1)) 
half <- max (t . sub*est . sub) 

while (half + 1 / N > 
N <- N + 2000 

b <- floor (sqrt(M)) # batch size 
t.sub <- qt(.9875, (N - b + 1)) 
est. sub <- subsainpling(chain[l :N] ) 
half <- max(t.sub*est.sub) 

} 
N 

half 

quant (chain[l:N]) 
t . sub*est . sub 
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3 T-Distribution Example 



Suppose our goal is to estimate the first two moments of a Students t distribution 
with 4 degrees of freedom and having density 



m{x) — 




Obviously, there is nothing about this that requires MCMC since we can easily 
calculate that EmX = and EmX^ = 2. Nevertheless, we will use a data 
augmentation algorithm based on the joint density 

SO that the full conditionals are X\y N(0, and Y\X ~ Gamma(5/2, 2 + 
/2). Consider the Gibbs sampler that updates X then Y so that a one-step 
transition looks like {x',y') — >■ {x,y') — >■ {x,y). 



3.1 Mcirkov Chain Sampler 

The following function samples from p iterations in the Markov chain starting 
from the value (1, 1). 

> t.gen <- functionCp, x = 1, y = 1) { 

+ loc <- length(x) 

+ junk <- double(p) 

+ X <- appendix, junk) 

+ 7 append (y, junk) 

+ for a in l:p) { 

+ j <- i + loc - 1 

+ x[(j + 1)1 <- rnormCl, 0, sqrt(l/y[j])) 

+ y[(j + 1)1 <- rgammad, 5/2, rate = (2 + x[(j + l)l~2/2)) 

+ } 

+ return (cbind(x, y)) 
+ } 



3.2 Calculaitons and Plots 

Using the sampler, we can run the simulation and create the necessary plots. 

> n <- 2000 

> iter <- seq(l, n) 

> set.seed(lOO) 

> sample <- t.gen((n - 1)) 

> X <- sample [, 1] 

> y <- sample [, 2] 

> RB <- 1/y 
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> x.mean <- cumsum(x) /seqCalong = x) 

> x2.mean <- cumsum(x~2)/seq(along = x) 

> RB.mean <- cumsum(l/y)/seq(along = y) 

> ts .plot (x. mean, main = "First Moment", xlab = "Iteration" , ylab = "", 
+ Iwd = 2, xlim = c(0, n)) 

> abline(h = 0) 

> abline(h = 0, Iwd = 2, Ity = 4) 

> legendC'topright", cC'Std", "RB"), Ity = c(l, 4), Iwd = 2) 
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> ts .plot Cx2. mean, main = "Second Moment", xlab = "Iteration", 
+ ylab = "", Iwd = 2, xlim = c(0, n)) 

> points Citer, ilB.mean, type = "1", Ity = 4, Iwd = 2) 

> ablineCh = 2) 

> legend("bottomright" , c("Std", "RB"), Ity = c(l, 4), Iwd = 2) 
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3.3 Interval Estimates 

The following chunk of code will calculate interval estimates using the same data 
from above. 

> obm.x <- sapplyCl :n, function(k) return(mcse(x[l :k] , meth = "OBM"))) 

> obm.RB <- sapplyCl :n, function(k) return(mcse(RB[l :k] , meth = "OBM"))) 

> obm.x2 <- sapplyCl :n, function returii(mcse(x[l:k]''2, meth = "OBM"))) 

Then we can plot the estimates for the first moment using standard methods. 
Notice, there is no uncertainty with the RB estimator in this setting. 

> upper <- x.mean + crit.obm * obm.x 

> lower <- x.mean - crit.obm * obm.x 

> ts. plot (x.mean, main = "First Moment", xlab = "Iteration" , ylab = "", 
+ Iwd = 2, xlim = c(0, n)) 

> ablineCh = 0) 

> abline(h = 0, Iwd = 2, Ity = 4) 

> legend("topright", c("Std", "RB"), Ity = c(l, 4), Iwd = 2) 

> points ("iter, upper, type = "1", Ity = 1, Iwd = 1) 

> points(iter, lower, type = "1", Ity = 1, Iwd = 1) 
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Iteration 

Now the plot for second moment for both the usual and RB estimator. 

> upper <- RB.mean + cr it. obm * obm.RB 

> lower <- RB.mean. - crit.obm * obm.RB 

> upper. x2 <- x2.meaii + crit.obm * obm.x2 

> lower. x2 <- x2.meaxi - crit.obm * obm.x2 

> ts.plot(x2.mesm, main = "Second Moment", xlab = "Iteration" , 

+ ylab = "", Iwd = 2, xlim = c(0, n) , ylim = c (min (lower. x2 [10 : n] ) , 
+ max (upper. x2[10:n]))) 

> points (iter , RB.mean, type = "1", Ity = 4, Iwd = 2) 

> abline(h = 2) 

> legend("bottomright" , c("Std", "RB") , Ity = c(l, 4), Iwd = 2) 

> points ("iter, upper, type = "I", Ity = 4, Iwd = 1) 

> points(iter, lower, type = "1", Ity = 4, Iwd = 1) 

> points (iter , upper. x2, type = "I", Ity = 1, Iwd = 1) 

> pointsdter, lower. x2, type = "1", Ity = 1, Iwd = 1) 
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Second Moment 




4 Estimating Marginals Example 

Suppose Yi\n, 9 ~ N(/x, 6) independently for i = 1, . . . , m where m > 3 and prior 
9) oc 9~^l'^. The target is the posterior density 

7r(M,e|t/) (X 0-('"+i)/2e-il(-'+te-M)^) 

where is the usual biased sample variance. It is easy to see that p\9^y ^ 
l:^{y,9/m) and that 9\fi,y ^ lG{{m-l)/2,m[s'^ + {y - iJ.)^]/2). We will conisder 
the Gibbs sampler that updates /x then 9 so that a one-step transition is given 
by in', 9') {f^,9') (Mi^)- We will use this Gibbs sampler to estimate the 
marginal densities of jj, and 9. 

4.1 Mcirkov Chain Sampler 

The first function provides an observation one step ahead using the Gibbs sam- 
pler. The second function results in p observations from the Markov chain. 

> ex2. gibbs <- function(m, t, y.bar =1, s2 = 4, n = 11) { 

+ alpha <- (n - l)/2 

+ beta <- n * (s2 + (y.bax - m)''2)/2 

+ t <- rgamma (1 , alpha, beta) 



19 



+ t <- 1/t 

+ m <- rnormCl , y.bar, sqrt (t/n) ) 

+ cbindCm, t) 
+ } 

> ex2.gen <- function (mc , p = 100, q = 1, y.bar = 1, s2 = 4, n = 11) { 

+ if (is .matrix (mc) == TRUE) { 

+ loc <- dim(mc) [1] 

+ } 

+ else { 

+ loc <- 1 

+ mc <- t (as. matrix (mc)) 

+ } 

+ junk <- matrix (rep (NA, p * 2) , ncol = 2) 

+ mc <- rbind(mc , junk) 

+ for (i in 1 :p) { 

+ j <- i + loc - 1 

+ mc[(j + 1) , 1 <- ex2 . gibbs (mc [j , 1] , mc[j, 2]) 

+ } 

+ return(mc) 
+ } 



4.2 Calculations and Settings 

Now suppose 171 = 11, y = 1 and = 4. We simulated 2000 realizations of the 
Gibbs sampler starting from (/ii, Ai) = (1, 1). The marginal density plots were 
created using the default settings for the density function in R. The bivariate 
density plot was created using R functions kde2d and persp. The resulting 
posterior is simple, so it is no surprise that the Gibbs sampler has been shown 



to converge in just a few iterations (Jones and Hobert, 2001 ). 



> n <- 2000 

> set.seed(lOO) 

> start <- c(l, 1) 

> sample <- ex2 . gen (start , (n - 1)) 

> mu <- sample [, 1] 

> theta <- sample [, 2] 

Now we can plot the marginal densities and corresponding histograms for /x 
and 9. 

> hist(mu, main = "", xlab = "", freq = F) 

> junk <- density (mu) 

> points (junk$x, junk$y, type = "1") 
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> hist(theta, main = "", xlab = "", freq = F, ylim = c(0, 0.18)) 

> junk <- density (theta) 

> points (junk$x, junk$y, type = "1") 



21 




Now we can plot the bivariate density plot. 

> library(MASS) 

> junk <- kde2d(mu, theta, n = 50, lims = c(-1.5, 3.5, 1, 15)) 

> perspCjunk, theta = 120, phi = 30, xlab = "mu" , ylab = "theta", 
+ zlab = "", expand =0.6) 
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4.3 Alternative Curve Estimation 



A clever technique for estimating a marginal is based on the same idea as RB 
estimators (Wei and Tanner 1990). The following chunk of code implements 



this for the current example for the marginal of /i. 

> hist(mu, main = "", xlab = "", freq = F) 

> junk <- density (mu) 

> points (junk$x, junk$y, type = "1") 

> avals <- seq(-3, 4, 0.01) 

> norm. avals <- dnorm(evals , 1, sqrt(mesm(theta)/ll)) 

> points (evals, norm. evals , type = "1", Ity = 4, Iwd = 2) 
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